R Notebook

Prepare the Workspace

knitr::opts_chunk$set(warning=FALSE, message=FALSE, error = FALSE) # rmd options
rm(list = ls()); invisible(gc()) # cleaning

Control Block

# f, f, 0.25 looks nice.
allowmigrants <- F # OPTIONS: T, F 
allowsympatry <- F # OPTIONS: T, F
minoverperc <- 0 # remove pairs that do not have thermal overlap (anagenesis)
costvar <- "ele"

Packages & Prefs

options(scipen = 999) # turn off scientific notation
'%notin%' <- Negate('%in%')
require(ggplot2) # load packages

Loading required package: ggplot2

require(GGally)

Loading required package: GGally
Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

require(viridis)

Loading required package: viridis
Loading required package: viridisLite

require(caper)

Loading required package: caper
Loading required package: ape
Loading required package: MASS
Loading required package: mvtnorm

library(dplyr)

Attaching package: ‘dplyr’

The following object is masked from ‘package:MASS’:

    select

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(stringr)
library(maps)
library(ape)
library(EnvStats)

Attaching package: ‘EnvStats’

The following object is masked from ‘package:MASS’:

    boxcox

The following objects are masked from ‘package:stats’:

    predict, predict.lm

The following object is masked from ‘package:base’:

    print.default

library(forecast)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

require(nlme)

Loading required package: nlme

Attaching package: ‘nlme’

The following object is masked from ‘package:forecast’:

    getResponse

The following object is masked from ‘package:dplyr’:

    collapse

require(geodist)

Loading required package: geodist

require(letsR)

Loading required package: letsR
Loading required package: raster
Loading required package: sp

Attaching package: ‘raster’

The following object is masked from ‘package:nlme’:

    getData

The following object is masked from ‘package:EnvStats’:

    cv

The following object is masked from ‘package:dplyr’:

    select

The following objects are masked from ‘package:MASS’:

    area, select

The following objects are masked from ‘package:ape’:

    rotate, zoom

require(spdep)

Loading required package: spdep
Loading required package: spData
To access larger datasets in this package, install the spDataLarge package with: `install.packages('spDataLarge', repos='https://nowosad.github.io/drat/',
type='source')`
Loading required package: sf
Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1
Registered S3 method overwritten by 'spdep':
  method   from
  plot.mst ape 

require(spatialreg)

Loading required package: spatialreg
Loading required package: Matrix
Registered S3 methods overwritten by 'spatialreg':
  method                   from 
  residuals.stsls          spdep
  deviance.stsls           spdep
  coef.stsls               spdep
  print.stsls              spdep
  summary.stsls            spdep
  print.summary.stsls      spdep
  residuals.gmsar          spdep
  deviance.gmsar           spdep
  coef.gmsar               spdep
  fitted.gmsar             spdep
  print.gmsar              spdep
  summary.gmsar            spdep
  print.summary.gmsar      spdep
  print.lagmess            spdep
  summary.lagmess          spdep
  print.summary.lagmess    spdep
  residuals.lagmess        spdep
  deviance.lagmess         spdep
  coef.lagmess             spdep
  fitted.lagmess           spdep
  logLik.lagmess           spdep
  fitted.SFResult          spdep
  print.SFResult           spdep
  fitted.ME_res            spdep
  print.ME_res             spdep
  print.lagImpact          spdep
  plot.lagImpact           spdep
  summary.lagImpact        spdep
  HPDinterval.lagImpact    spdep
  print.summary.lagImpact  spdep
  print.sarlm              spdep
  summary.sarlm            spdep
  residuals.sarlm          spdep
  deviance.sarlm           spdep
  coef.sarlm               spdep
  vcov.sarlm               spdep
  fitted.sarlm             spdep
  logLik.sarlm             spdep
  anova.sarlm              spdep
  predict.sarlm            spdep
  print.summary.sarlm      spdep
  print.sarlm.pred         spdep
  as.data.frame.sarlm.pred spdep
  residuals.spautolm       spdep
  deviance.spautolm        spdep
  coef.spautolm            spdep
  fitted.spautolm          spdep
  print.spautolm           spdep
  summary.spautolm         spdep
  logLik.spautolm          spdep
  print.summary.spautolm   spdep
  print.WXImpact           spdep
  summary.WXImpact         spdep
  print.summary.WXImpact   spdep
  predict.SLX              spdep

Attaching package: ‘spatialreg’

The following objects are masked from ‘package:spdep’:

    anova.sarlm, as_dgRMatrix_listw, as_dsCMatrix_I, as_dsCMatrix_IrW, as_dsTMatrix_listw, as.spam.listw, bptest.sarlm, can.be.simmed, cheb_setup, coef.gmsar,
    coef.sarlm, coef.spautolm, coef.stsls, create_WX, deviance.gmsar, deviance.sarlm, deviance.spautolm, deviance.stsls, do_ldet, eigen_pre_setup, eigen_setup,
    eigenw, errorsarlm, fitted.gmsar, fitted.ME_res, fitted.sarlm, fitted.SFResult, fitted.spautolm, get.ClusterOption, get.coresOption, get.mcOption,
    get.VerboseOption, get.ZeroPolicyOption, GMargminImage, GMerrorsar, griffith_sone, gstsls, Hausman.test, HPDinterval.lagImpact, impacts, intImpacts, Jacobian_W,
    jacobianSetup, l_max, lagmess, lagsarlm, lextrB, lextrS, lextrW, lmSLX, logLik.sarlm, logLik.spautolm, LR.sarlm, LR1.sarlm, LR1.spautolm, LU_prepermutate_setup,
    LU_setup, Matrix_J_setup, Matrix_setup, mcdet_setup, MCMCsamp, ME, mom_calc, mom_calc_int2, moments_setup, powerWeights, predict.sarlm, predict.SLX, print.gmsar,
    print.ME_res, print.sarlm, print.sarlm.pred, print.SFResult, print.spautolm, print.stsls, print.summary.gmsar, print.summary.sarlm, print.summary.spautolm,
    print.summary.stsls, residuals.gmsar, residuals.sarlm, residuals.spautolm, residuals.stsls, sacsarlm, SE_classic_setup, SE_interp_setup, SE_whichMin_setup,
    set.ClusterOption, set.coresOption, set.mcOption, set.VerboseOption, set.ZeroPolicyOption, similar.listw, spam_setup, spam_update_setup, SpatialFiltering,
    spautolm, spBreg_err, spBreg_lag, spBreg_sac, stsls, subgraph_eigenw, summary.gmsar, summary.sarlm, summary.spautolm, summary.stsls, trW, vcov.sarlm, Wald1.sarlm

require(rnaturalearth)

Loading required package: rnaturalearth

require(rnaturalearthdata)

Loading required package: rnaturalearthdata

require(rgeos)

Loading required package: rgeos
rgeos version: 0.5-5, (SVN revision 640)
 GEOS runtime version: 3.8.1-CAPI-1.13.3 
 Linking to sp version: 1.4-2 
 Polygon checking: TRUE 

require(sf)
require(rgdal)

Loading required package: rgdal
rgdal: version: 1.5-18, (SVN revision 1082)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.1.4, released 2020/10/20
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/sf/gdal
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/proj
Linking to sp version:1.4-4
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

require(raster)
world <- ne_coastline(scale = "medium", returnclass = "sf")
vlog <- function(x){
   log( x + abs(min( x , na.rm = T)) + 1)
}

setwd("~/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Analyze_Processed_Cluster_Outputs/Data")

The working directory was changed to /Users/boterolab1/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Analyze_Processed_Cluster_Outputs/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

exclusion <- read.csv(file = "exclusion_nonsimpatric_nonmigrant.csv")
exclusion$realm1red[is.na(exclusion$realm1red)] <- "NA" # NA is north america not R's NA value. fix.
exclusion$realm2green[is.na(exclusion$realm2green)] <- "NA"

Load Main Data

# main dataframe ---------------------------------
setwd("~/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Process_Cluster_Outputs/Data")
load(file =  "Pair_Barrier_Data_FEB2021.rdata")
mydata <- mypairdata; rm(mypairdata)
rownames(mydata) <- mydata$Species.1
mydatahold <- mydata

initial masks

# migration ---
if(allowmigrants == F){
  mydata <- mydata[which(mydata$Migration == 1.0),]
}

# patry
if(allowsympatry == F){
  mydata <- mydata[which(mydata[,paste0(costvar, "_c0")] > 0 ),] # doesnt matter if you use ele, mat, vart paths here, will be same answer
}

sort

mydata$uniquePairId == exclusion$mydata.uniquePairId

  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [31] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [91] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[151] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[211] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[241] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

x <- mydata[which(mydata$Species.1 == "Apteryx_owenii"),]
mydata <- mydata[order(match(mydata$uniquePairId,exclusion$mydata.uniquePairId)), ]
y <- mydata[which(mydata$Species.1 == "Apteryx_owenii"),]
# sum(y!=x, na.rm = T) # checks.
# head(mydata)
mydata$uniquePairId == exclusion$mydata.uniquePairId

  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [37] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [73] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[109] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[145] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[181] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[217] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[253] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

rm(x,y)

# # Basic range maps for all pairs (no paths)
# wdPAM <- "/Users/boterolab1/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/PREP/PAM/Data"
# setwd(wdPAM); load("cbPAM.rdata")
# setwd("~/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Process_Cluster_Outputs/Data")
# pdf("pairmaps2.pdf", width = 19, height = 9.25)
# for (i in 1:nrow(mydata)){
# x <- cbPAM[,c("Longitude(x)","Latitude(y)",mydata$Species.1bl[i])]
# x <- as.data.frame(x[x[,3] == 1,])
# if(ncol(x) == 1) {
#   x <- t(x)
#   colnames(x) <- c("lon", "lat", "pres")
#   x <- as.data.frame(x)
# } else {
#   colnames(x) <- c("lon", "lat", "pres")
# }
# 
# y <- cbPAM[,c("Longitude(x)","Latitude(y)",mydata$Species.2bl[i])]
# y <- as.data.frame(y[y[,3] == 1,])
# if(ncol(y) == 1) {
#   y <- t(y)
#   colnames(y) <- c("lon", "lat", "pres")
#   y <- as.data.frame(y)
# 
# } else {
#   colnames(y) <- c("lon", "lat", "pres")
# }
# z <- ggplot(world)+
#     geom_sf() +
#     geom_point(data = x, aes(y=lat, x=lon), color = "red") +
#     geom_point(data = y, aes(y=lat, x=lon), color = "green") +
#     theme_bw() +
#     ggtitle(i)
# print(z)
# }
# print(i)
# dev.off()

mydata_exclusion <- cbind(mydata, exclusion)
save(mydata_exclusion, file = "mydata_exclusion.rdata")
mydata$realm1 <- exclusion$realm1red
mydata$realm2 <- exclusion$realm2green
mydata$realm <- paste0(mydata$realm1, mydata$realm2)
table(mydata$realm)

AAAA AAIM ATAT ATIM ATPA IMAA IMAT IMIM IMNT NANA NTAA NTNA NTNT OCOC PAAT PANA PAPA 
  30    1   45    3    1    3    3   27    1    4    1    1  135    2    3    2    6 

mydata$realm[mydata$realm == "AAIM"]; mydata$realm[mydata$realm == "IMAA"] <-  "AAIM"

[1] "AAIM"

mydata$realm[mydata$realm == "ATIM"]; mydata$realm[mydata$realm == "IMAT"] <-  "ATIM"

[1] "ATIM" "ATIM" "ATIM"

mydata$realm[mydata$realm == "ATPA"]; mydata$realm[mydata$realm == "PAAT"] <-  "ATPA"

[1] "ATPA"

mydata$realm[mydata$realm == "IMNT"]; mydata$realm[mydata$realm == "NTIM"] <-  "IMNT"

[1] "IMNT"

mydata$realm[mydata$realm == "NTAA"]; mydata$realm[mydata$realm == "AANT"] <-  "NTAA"

[1] "NTAA"

mydata$realm[mydata$realm == "NTNA"]; mydata$realm[mydata$realm == "NANT"] <-  "NTNA"

[1] "NTNA"

mydata$realm[mydata$realm == "PANA"]; mydata$realm[mydata$realm == "NAPA"] <-  "PANA"

[1] "PANA" "PANA"

mydata$realm <- as.factor(mydata$realm); mydata$realm <- relevel(mydata$realm, "NTNT")

mydata$landgap <- as.logical(exclusion$island)

mydata$cosmopolitan <- as.logical(exclusion$cosmopolitan)

mydata$new.old <- as.logical(exclusion$new.old)

rm(exclusion, mydata_exclusion)

Impose masks & do calcuations

# filter cosmopolitan and new/old world species (there are relatively few after imposing previous masks.)
mydata <- mydata[which(mydata$cosmopolitan == FALSE & mydata$new.old == FALSE),]

# dependent variable: elevational barrier size ---
mydata$cost <- mydata[, paste0(costvar, "_c25")] 

# data filtering -----------------------
#  thermal overlap ---
mydata$MAT_overlap <- mydata[,paste0("MAT", "_ov_perc_smrnge")]
mydata <- mydata[mydata$MAT_overlap > minoverperc,]

# update sort order --------------------
mydata$sortorder <- seq(1:nrow(mydata))

# longitude ----------------------------
mydata$lon <- mydata[,paste0("lon_mean_pair_", costvar, "_c25")]

# latitude ----------------------------
mydata$lat <- mydata[,paste0("lat_mean_pair_", costvar, "_c25")]

# temperature breadth -----------------
mydata$tas_breadth <- mydata$tas_range # mean(mean(sp1 annual tas range -- one value per cell), mean(sp2 annual tas range -- one value per cell))

# mean annual temperature --------------
mydata$tas_position <- mydata$tas_mean # mean(mean(sp1 annual tas mean  -- one value per cell), mean(sp2 annual tas mean  -- one value per cell))

# precipitation breadth --------------- 
mydata$pcp_breadth <- mydata$pcp_range # mean(mean(sp1 annual pcp range  -- one value per cell), mean(sp2 annual pcp range  -- one value per cell))

# precipitation breadth --------------- 
mydata$pcp_position <- mydata$pcp_mean # mean(mean(sp1 annual pcp mean  -- one value per cell), mean(sp2 annual pcp mean  -- one value per cell))

# distance -----------------------------
mydata$distance <- mydata[,paste0("centroid_distance_",costvar,"_c25")]

# mountain mass ------------------------
mtns <- readOGR(dsn="~/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Other_Input_Data/GMBA", layer="GMBA Mountain Inventory_v1.2-World", verbose = FALSE)
wdPAM <- "/Users/boterolab1/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/PREP/PAM/Data"
setwd(wdPAM); load("LonLat_BirdPAM_raster.rdata")

The working directory was changed to /Users/boterolab1/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/PREP/PAM/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

mtns <- rasterize(mtns, LonLat_BirdPAM_raster)
mtns@data@values[!is.na(mtns@data@values)] <- 1 # replace mountain IDs with simple coding. 1 for mountain...
mtns@data@values[is.na(mtns@data@values)] <- 0 # ...0 for no mountain

mydata$mtn_mass <- NA
for (i in 1:nrow(mydata)) {
  coords1<-data.frame(lon=mydata$lon[i], lat=mydata$lat[i]); coordinates(coords1)<-c("lon","lat"); crs(coords1)<-crs(LonLat_BirdPAM_raster) # get coordinates
  z <- extract(mtns, coords1, buffer = raster::pointDistance(c(0,0), c(0,8), lonlat = T))
  mydata$mtn_mass[i] <- sum(z[[1]]) / length(z[[1]])
}

raster::pointDistance(c(0,0), c(0,8), lonlat = T) # corresponds to a radius of just about 8 degrees.

[1] 884651.7

raster::pointDistance(c(75,75), c(75,(75+8)), lonlat = T) # polar circles are bigger, but not that much bigger. so should be OK.

[1] 893209.5

raster::plot(mtns)
raster::plot(world$geometry, add = T)

plotrix::draw.circle(70, 40, 8, nv = 1000, border = "hotpink", lty = 1, lwd = 1)
plotrix::draw.circle(140, -5, 8, nv = 1000, border = "hotpink", lty = 1, lwd = 1)

plotrix::draw.circle(-80, 10, 8, nv = 1000, border = "hotpink", lty = 1, lwd = 1)
plotrix::draw.circle(-70, -50, 8, nv = 1000, border = "hotpink",lty = 1, lwd = 1)

x <- ggplot(world)+
  geom_sf() +
  geom_point(data = mydata[order(mydata[, "mtn_mass"], decreasing = F),], aes(y=lat, x=lon, color = mtn_mass), alpha = 0.9) +
  ggtitle(i)+scale_color_viridis()
print(x)

# water buffering ----------------------
wtr <- read_sf("~/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Other_Input_Data/ne_50m_ocean/ne_50m_ocean.shp")
temp <- mtns
values(temp) <- 100
temp <- mask(temp, wtr)
wtr <- mtns
values(wtr) <- 0
wtr[temp == 100] <- 1; rm(temp)
plot(wtr)

mydata$wtrcells <- NA
mydata$wtrmass <- NA
for (i in 1:nrow(mydata)) {
  coords1<-data.frame(lon=mydata$lon[i], lat=mydata$lat[i]); coordinates(coords1)<-c("lon","lat"); crs(coords1)<-crs(LonLat_BirdPAM_raster) # get coordinates
  z <- extract(wtr, coords1, buffer = raster::pointDistance(c(0,0), c(0,8), lonlat = T))
  mydata$miniplot[i] <- length(z[[1]])
  mydata$wtrcells[i] <- sum(z[[1]])
  mydata$wtrmass[i] <- sum(z[[1]]) / length(z[[1]])
}
mydata$water_buffering <- mydata$wtrmass

x <- ggplot(world)+
  geom_sf() +
  geom_point(data = mydata[order(mydata[, "water_buffering"], decreasing = F),], aes(y=lat, x=lon, color = water_buffering), alpha = 0.9) +
  ggtitle(i)+scale_color_viridis()
print(x)

# dispersal ability --------------------
dispab <- read.csv("~/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Other_Input_Data/Bird Hand-Wing Index/Dataset HWI 2020-04-10.csv")
mydata$dispersal_ability <- NA
for (i in 1:nrow(mydata)){
  dispab_sp1 <- dispab$HWI[dispab$IUCN.name == mydata$Species.1bl[i]]
  dispab_sp2 <- dispab$HWI[dispab$IUCN.name == mydata$Species.2bl[i]]
  dispab_pair <- mean(c(dispab_sp1, dispab_sp2), na.rm = T)
  mydata$dispersal_ability[i] <- dispab_pair
  rm(dispab_sp1, dispab_sp2, dispab_pair)
}
mydata <- mydata[!is.nan(mydata$dispersal_ability),]

# pair age -----------------------------
mydata$pair_age <- mydata$Pair.age..MY.

# length of boundary -------------------
mydata$boundary_length <- mydata$boundary_length_ele_c25

# retain cols of interest only.
mydata <- mydata[,c("uniquePairId", "Species.1", "Species.2", "Species.1bl", "Species.2bl", "cost", "lat", "lon",
                    "tas_breadth","tas_position", "pcp_breadth","pcp_position", "mtn_mass", "water_buffering", 
                    "dispersal_ability", "pair_age", "distance", "boundary_length", "MAT_overlap", "realm",
                    "landgap")]

load phylo and prune

sum(mydata$Species.1 != tree$tip.label) # sorted (but still specify form below for safety)

[1] 0

Save freame for CB

setwd("~/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Analyze_Processed_Cluster_Outputs/Data")

The working directory was changed to /Users/boterolab1/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Analyze_Processed_Cluster_Outputs/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

save(mydata, file = "ele_data_for_CB.rdata")
write.csv(mydata, file = "ele_data_for_CB.csv")

Neighbors for sptial analysis.

distm <- mydata[,c("lon", "lat")]
distm <- geodist(distm, measure = "geodesic")
rownames(distm) <- rownames(mydata)
colnames(distm) <- rownames(mydata)
basecols <- ncol(mydata)

# neightbors
coords<-cbind(mydata$lon, mydata$lat); coords<-as.matrix(coords) ; row.names(coords)<-rownames(mydata)
k1 <- knn2nb(knearneigh(coords, longlat = T))

knearneigh: identical points found

nb<- dnearneigh(coords,row.names = row.names(coords), d1=0,d2=max(unlist(nbdists(k1, coords, longlat = T))),longlat=T)

plots for transformations

for (i in c("cost", "lat", "lon","tas_breadth","tas_position","pcp_breadth","pcp_position", "mtn_mass", "dispersal_ability", "pair_age", "distance", "boundary_length", "MAT_overlap")) {
hist(mydata[, i], breaks = 50, main = i)
}

world plots

setwd("~/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Analyze_Processed_Cluster_Outputs/Data")

The working directory was changed to /Users/boterolab1/Box Sync/CB_VF_Shared/Dry_Lab/Projects/JMPH/Analyze_Processed_Cluster_Outputs/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

pdf(file="variable_maps.pdf", width = 10, height = 7)
for (i in c("cost", "lat", "lon","tas_breadth","tas_position","pcp_breadth","pcp_position", "mtn_mass",  "water_buffering", "dispersal_ability", "pair_age", "distance", "boundary_length", "MAT_overlap", "realm", "landgap")) {
  if(i %in% c("cost")){
    myc <- mydata[, i]; myc <- vlog(myc)
    x <- ggplot(world)+
      geom_sf() +
      geom_point(data = mydata[order(mydata[, i], decreasing = F),], aes(y=lat, x=lon, color = myc), alpha = 0.85) +
      ggtitle(paste0("ln ",i))+
      scale_color_viridis()
    print(x)
  } else if (i %in% c("realm", "landgap")){
    x <- ggplot(world)+
      geom_sf() + 
      geom_point(data = mydata[order(mydata[, "cost"], decreasing = F),], aes(y=lat, x=lon, color = get(i)), alpha = 0.85) +
      ggtitle(paste0(i))+
      scale_color_viridis(discrete = T)
    print(x) 
  } else {
    x <- ggplot(world)+
      geom_sf() + 
      geom_point(data = mydata[order(mydata[, "cost"], decreasing = F),], aes(y=lat, x=lon, color = get(i)), alpha = 0.85) +
      ggtitle(i)+
      scale_color_viridis()
    print(x) 
  }
}
dev.off()

null device 
          1 

# vector 1 added ---
m <- gls(scale(I(vlog(cost))) ~ scale(I(vlog(tas_breadth))) + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30, correlation = corPagel(0.99, phy = tree,fixed = F, form = ~Species.1), data = mydata, method = "REML"); summary(m); cor(predict(m),scale(I(vlog(mydata$cost))))

Generalized least squares fit by REML
  Model: scale(I(vlog(cost))) ~ scale(I(vlog(tas_breadth))) + V22 + V23 +      V24 + V25 + V26 + V27 + V28 + V29 + V30 
  Data: mydata 

Correlation Structure: corPagel
 Formula: ~Species.1 
 Parameter estimate(s):
   lambda 
0.1132853 

Coefficients:

 Correlation: 
                            (Intr) s(I((_ V22    V23    V24    V25    V26    V27    V28    V29   
scale(I(vlog(tas_breadth))) -0.030                                                               
V22                          0.060 -0.016                                                        
V23                         -0.033  0.002 -0.008                                                 
V24                         -0.031  0.016 -0.012  0.005                                          
V25                          0.015  0.020 -0.003 -0.014  0.004                                   
V26                          0.032  0.007  0.004 -0.015 -0.006  0.006                            
V27                         -0.018 -0.020  0.003  0.003  0.001 -0.010 -0.003                     
V28                          0.068  0.000  0.001  0.000  0.002  0.003 -0.017  0.000              
V29                         -0.038  0.001 -0.012  0.002  0.017 -0.001 -0.018  0.001  0.000       
V30                          0.031  0.024 -0.007 -0.013 -0.011  0.013  0.061 -0.015  0.001 -0.014

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-3.0866866 -0.5483974  0.1080253  0.6419830  2.6354891 

Residual standard error: 0.8999846 
Degrees of freedom: 235 total; 224 residual
          [,1]
[1,] 0.5013446

matx <- as.matrix(m$residuals); rownames(matx) <- rownames(mydata)
spac <- lets.correl(x=matx, y=distm, z=12, equidistant = T, plot = T)

moran.test(residuals(m), nb2listw(nb))$p.value # too much.

[1] 0.9215695

hist(resid(m))

plot(m, resid(., type = "p") ~ fitted(.), abline = 0)

plot(m, scale(I(vlog(cost))) ~ fitted(.), abline = c(0,1))

qqnorm(m)

for(i in c("V22", "V23", "V24", "V25", "V26", "V27", "V28", "V29", "V30")){
  x <- ggplot(world)+
  geom_sf() +
  geom_point(data = mydata[order(mydata[,i]),], aes(y=lat, x=lon, color = mydata[,i]), alpha = 0.9) +
  scale_color_viridis()+
  ggtitle(i)
  print(x)
}

Is thermal niche breadth predicted by latitude plus latitude2?

# vector 1 added ---
m <- gls(scale(I(vlog(tas_breadth))) ~ scale(lat) + scale(I(lat^2)) + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30 + V31 + V32 + V33 + V34 + V35 + V36 + V37, correlation = corPagel(0.99, phy = tree,fixed = F, form = ~Species.1), data = mydata, method = "REML"); summary(m); cor(predict(m),scale(I(vlog(mydata$cost))))

Generalized least squares fit by REML
  Model: scale(I(vlog(tas_breadth))) ~ scale(lat) + scale(I(lat^2)) +      V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30 + V31 +      V32 + V33 + V34 + V35 + V36 + V37 
  Data: mydata 

Correlation Structure: corPagel
 Formula: ~Species.1 
 Parameter estimate(s):
   lambda 
0.5675592 

Coefficients:

 Correlation: 
                (Intr) scl(l) s(I(^2 V22    V23    V24    V25    V26    V27    V28    V29    V30    V31    V32    V33    V34    V35    V36   
scale(lat)       0.045                                                                                                                       
scale(I(lat^2)) -0.077  0.055                                                                                                                
V22              0.011 -0.056 -0.055                                                                                                         
V23              0.070  0.021  0.019  0.048                                                                                                  
V24              0.053  0.032 -0.058  0.076  0.005                                                                                           
V25             -0.024  0.005 -0.134  0.208  0.006  0.167                                                                                    
V26             -0.072 -0.037  0.034  0.040 -0.044 -0.007 -0.014                                                                             
V27             -0.007 -0.004 -0.078  0.050 -0.032 -0.022  0.090  0.028                                                                      
V28              0.028 -0.005 -0.006 -0.069 -0.003  0.037 -0.018 -0.057 -0.047                                                               
V29              0.002  0.011 -0.022 -0.006  0.002 -0.052  0.005  0.011  0.048 -0.004                                                        
V30             -0.038 -0.018 -0.018 -0.017 -0.009 -0.013 -0.005  0.030 -0.011 -0.008 -0.002                                                 
V31              0.005 -0.030 -0.009 -0.062 -0.030  0.021 -0.006 -0.046 -0.019  0.014 -0.023 -0.004                                          
V32              0.001 -0.005  0.012  0.022  0.004 -0.002 -0.021  0.004 -0.003 -0.002 -0.008 -0.011 -0.037                                   
V33             -0.008 -0.016 -0.070  0.043 -0.012 -0.014  0.026  0.012  0.050 -0.015  0.014  0.003 -0.007 -0.003                            
V34             -0.021  0.062  0.018 -0.018  0.010 -0.030 -0.030  0.026  0.022 -0.022  0.074  0.012 -0.027  0.010  0.036                     
V35             -0.005  0.012 -0.007 -0.006 -0.010 -0.013 -0.012 -0.001 -0.007 -0.001 -0.001  0.001 -0.020  0.002  0.003  0.004              
V36             -0.045 -0.014 -0.034  0.063 -0.007 -0.005  0.055  0.025  0.041 -0.022  0.021  0.008 -0.073 -0.028  0.013 -0.004 -0.003       
V37              0.027 -0.017 -0.021  0.004  0.017  0.013 -0.006  0.010 -0.006 -0.003 -0.015  0.002  0.008 -0.014 -0.005 -0.012  0.001 -0.008

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-2.0731629 -0.6096382 -0.1790288  0.2623287  3.2865813 

Residual standard error: 0.5443762 
Degrees of freedom: 235 total; 216 residual
          [,1]
[1,] 0.1123255

matx <- as.matrix(m$residuals); rownames(matx) <- rownames(mydata)
spac <- lets.correl(x=matx, y=distm, z=12, equidistant = T, plot = T)

moran.test(residuals(m), nb2listw(nb))$p.value # too much.

[1] 0.9484634

hist(resid(m))

plot(m, resid(., type = "p") ~ fitted(.), abline = 0)

plot(m, scale(I(vlog(cost))) ~ fitted(.), abline = c(0,1))

qqnorm(m)

for(i in c("V22", "V23", "V24", "V25", "V26", "V27", "V28", "V29", "V30", "V31", "V32", "V33", "V34", "V35", "V36", "V37")){
  x <- ggplot(world)+
  geom_sf() +
  geom_point(data = mydata[order(mydata[,i]),], aes(y=lat, x=lon, color = mydata[,i]), alpha = 0.9) +
  scale_color_viridis()
  print(x)
}

Is cost predicted by latitude plus latitude2?

# vector 1 added ---
m <- gls(scale(I(vlog(cost))) ~ scale(lat) + scale(I(lat^2)) + V22 + V23 + V24 + V25, correlation = corPagel(0.99, phy = tree,fixed = F, form = ~Species.1), data = mydata, method = "REML"); summary(m); cor(predict(m),scale(I(vlog(mydata$cost))))

Generalized least squares fit by REML
  Model: scale(I(vlog(cost))) ~ scale(lat) + scale(I(lat^2)) + V22 + V23 +      V24 + V25 
  Data: mydata 

Correlation Structure: corPagel
 Formula: ~Species.1 
 Parameter estimate(s):
   lambda 
0.1048059 

Coefficients:

 Correlation: 
                (Intr) scl(l) s(I(^2 V22    V23    V24   
scale(lat)       0.001                                   
scale(I(lat^2)) -0.064  0.104                            
V22             -0.070 -0.017  0.013                     
V23             -0.041 -0.007 -0.006  0.011              
V24             -0.031  0.013  0.018  0.009  0.001       
V25              0.032  0.005 -0.013 -0.021 -0.008 -0.004

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-3.2135868 -0.5833627  0.1685737  0.6876315  2.6281318 

Residual standard error: 0.9181122 
Degrees of freedom: 235 total; 228 residual
          [,1]
[1,] 0.4518531

matx <- as.matrix(m$residuals); rownames(matx) <- rownames(mydata)
spac <- lets.correl(x=matx, y=distm, z=12, equidistant = T, plot = T)

moran.test(residuals(m), nb2listw(nb))$p.value # too much.

[1] 0.7738208

hist(resid(m))

plot(m, resid(., type = "p") ~ fitted(.), abline = 0)

plot(m, scale(I(vlog(cost))) ~ fitted(.), abline = c(0,1))

qqnorm(m)

for(i in c("V22", "V23", "V24", "V25")){
  x <- ggplot(world)+
  geom_sf() +
  geom_point(data = mydata[order(mydata[,i]),], aes(y=lat, x=lon, color = mydata[,i]), alpha = 0.9) +
  scale_color_viridis()
  print(x)
}

Sensitivity Analyses

# 1 Pair age (all v. < 8mya (end of uplift of Andes))
# 2 Distance (all v. < 1500*1000) (1500 / 110 = ~ 22 degrees)
# 3 MAT_overlap (> 0% v. > 75% (more restrictive == more conservative for this measure.))
# 4 landgap (all v. nogap) *ALL GAPS ARE < 110km (two water grid cells marked as land for having >50% land @ 0.5 degree resolution.)